arXiv: 1503.0388lv2 [physics.plasm-ph] 16 Mar 2015 


Bump formation in the runaway electron tail 

J. Decker/’* E. Hirvijoki/ O. Embreus/ Y. Peysson/ A. Stahl/ 1. Pusztai/ and T. Fiilop^ 

^Ecole Polytechnique Federate de Lausanne (EPFL), Centre de Recherches 
en Physique des Plasmas (CRPP), CH-1015 Lausanne, Switzerland 
^Department of Applied Physics, Chalmers University of Technology, SE-41296 Gothenburg, Sweden 
^ CEA, IRFM, F-13108 Saint-Paul-lez-Durance, Franee 
(Dated: March 17, 2015) 

Runaway electrons are generated in a magnetized plasma when the parallel electric field exceeds 
a critical value. For such electrons with energies typically reaching tens of MeV, the Abraham- 
Lorentz-Dirac (ALD) radiation force, in reaction to the synchrotron emission, is significant and can 
be the dominant process limiting the electron acceleration. The effect of the ALD-force on runaway 
electron dynamics in a homogeneous plasma is investigated using the relativistic finite-difference 
Fokker-Planck codes LUKE [Decker & Peysson, Report EUR-CEA-FC-1736, Euratom-CEA, (2004)] 
and CODE [Landreman et al, Comp. Phys. Comm. 185, 847 (2014)]. Under the action of the ALD 
force, we find that a bump is formed in the tail of the electron distribution function if the electric field 
is sufficiently large. We also observe that the energy of runaway electrons in the bump increases with 
the electric field amplitude, while the population increases with the bulk electron temperature. The 
presence of the bump divides the electron distribution into a runaway beam and a bulk population. 
This mechanism may give rise to beam-plasma types of instabilities that could in turn pump energy 
from runaway electrons and alter their confinement. 


I. INTRODUCTION 


Runaway electrons are typically generated in plasmas 
in the presence of large electric fields E > Ec, where the 
critical field Ec is defined as [1] 


Ec 
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where n is the electron density, m is the electron rest 
mass, c is the speed of light, e is the elementary charge, 
and In A is the Coulomb logarithm. 

In connection with the sudden cooling in a tokamak 
disruption, a strong electric field is induced, which leads 
to the generation of a large number of runaway electrons. 
In certain cases a significant fraction of the initial toroidal 
current can be driven by a beam of runaway electrons. 
The formation of such energetic runaway beams would 
represent a serious threat for reactor-size machines such 
as ITER]2]. Consequently, a considerable research ef¬ 
fort is currently undertaken to prevent the formation of 
large runaway beams during tokamak disruptions, or to 
design a controlled damping scenario for runaway beams 
if they cannot be avoided [2-6]. The condition E > Ec 
for runaway electron generation can also be met during 
the plasma start-up or ramp-down. During the flat-top 
phase, runaways can appear if the density is sufficiently 
low in Ohmic plasmas (since Ec oc n), or if an externally 
applied source of current is suddenly modified. 

Experimental measurements show that the maximum 
runaway electron energy does not increase indefinitely 
with time but instead reaches a limit in the tens of MeV 
range [7]. One of the possible mechanisms that could 
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provide an explanation for this limit is the Abraham- 
Lorentz-Dirac (ALD) radiation force [8] in reaction to 
the synchrotron emission due to the particle motion in 
a magnetic field. For electrons in the MeV range, this 
force can be significant and contribute to limit the par¬ 
ticle acceleration [9]. The synchrotron emission, which 
carries energy away from the electrons, is also used as 
a diagnostic tool for the runaway population [10]. The 
ALD force is characterized by the synchrotron radiation 
reaction time scale Tr, given by 


67r£o (mc)^’ 

where B is the magnetic field. 

In the present paper, the runaway electron dynam¬ 
ics in a homogeneous plasma is investigated using the 
relativistic finite-difference guiding-center Fokker-Planck 
codes LUKE [11, 12] and CODE [13, 14]. The elec¬ 
tron distribution function evolves under the combined 
influence of Coulomb collisions, electric field accelera¬ 
tion, and the ALD radiation reaction force. Under a con¬ 
stant parallel electric field (with respect to the magnetic 
field) A|| > Ec, the electron distribution never reaches 
a steady-state in the absence of ALD radiation reaction 
force. Conversely, it is shown in Sec. Ill that when the 
effect of the ALD force is included, the electron distribu¬ 
tion evolves towards a steady-state solution. This solu¬ 
tion exists even though the synchrotron emission vanishes 
for electrons with purely parallel motion in a uniform 
magnetic field. In fact, the expansion of the electron dis¬ 
tribution towards higher energies is limited by collisional 
pitch-angle scattering, which is enhanced by the strong 
perpendicular anisotropy arising from the combination of 
electric field acceleration and synchrotron radiation re¬ 
action force. This process is found to limit the runaway 
electron population to energies far below the value for 
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which the contribution from the magnetic field curvature 
to the ALD radiation reaction force becomes significant 
[9]. Consequently, it is justified to use the homogeneous 
plasma limit to study the dynamics of runaway electrons 
in the core region of tokamaks. 

In addition, we show that if the electric field ampli¬ 
tude is sufficiently large, a bump appears in the runaway 
electron tail of the steady-state distribution function, in 
accordance with analytical predictions [15]. This bump, 
which peaks on the parallel axis in momentum space, is 
entirely located in the runaway region. The steady-state 
population of electrons in the bump is found to increase 
with the bulk electron temperature Tg, while their aver¬ 
age energy increases with the electric field amplitude 
and decreases with the amplitude of the ALD radiation 
reaction force, which is proportional to B^. For certain 
parameters, the bump in the electron distribution tail 
encompasses almost the entire runaway electron popula¬ 
tion, thus formally dividing the distribution into a bulk 
population and a runaway beam. 

The implementation of the synchrotron reaction force 
in the kinetic equation is described in Sec. II. The time 
evolution of the electron distribution calculated by the 
Fokker-Planck modelling code LUKE is presented in 
Sec. III. The properties of the steady-state distribution 
function and the mechanism leading to the formation of 
a bump are described in Sec. IV, where a comparison be¬ 
tween the codes LUKE and CODE is also presented. The 
bump is characterized in Sec. V as a function of the elec¬ 
tric field amplitude, ALD radiation reaction, bulk elec¬ 
tron temperature, and ion effective charge. Implications 
of the ALD radiation reaction force and bump-in-tail for¬ 
mation are discussed in the Conclusions, Sec. VI. 


II. SYNCHROTRON REACTION FORCE IN 
THE KINETIC EQUATION 

A. Kinetic equation for charged particles in a 
magnetized plasma 


mass m is given by 


Bf B B 


cle species a and b (including intra-species collisions) 
and (x, p) are the equations of motion associated with 
phase-space coordinates (x, p). Here x is the particle 
position and p = ymv is the particle momentum, with 
7 = 1/ -^/l — v'^/c^ = \/l + I (rncY the relativistic fac¬ 
tor. In the Fokker-Planck limit, the Coulomb collision 
operator is given by 


CMfa.h] = - 
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where 'K.abifb] is the collisional friction vector and 
'Oabifb] is the collisional diffusion tensor. The relativistic 
Braams-Karney collision operator is used in this paper 
[16, 17]. 

So-called knock-on collisions represent a 1/ In A correc¬ 
tion to the collision operators. However, when the run¬ 
away population becomes significant, these collisions can 
play an important role as they give rise to an avalanche 
effect that can significantly increase the runaway growth 
rate. This secondary runaway generation is neglected in 
the present work, which is restricted to situations where 
the runaway population is sufficiently small for secondary 
electron generation to be negligible. However, it is pos¬ 
sible that the ALD radiation reaction force has a signif¬ 
icant effect on the secondary runaway generation. Such 
considerations will be the subject of future work. 

The equations of motion combine the Hamiltonian mo¬ 
tion from the electric and magnetic fields E and B, and 
the effect of the ALD radiation reaction Fald : 

X =v, (5) 

P =9 (E -I- V X B) -I- Fald = F^; -I- F^ + Fald- (6) 

The Abraham-Lorentz-Dirac force describes momentum 
loss in reaction to the synchrotron radiation, and takes 
the form [8] 


Fald = 
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In magnetically confined fusion plasmas, the magnetic 
force Fm = gv x B characterized by the Larmor frequency 
ujc = qB/m typically dominates both the electric force 
Fb = gE and the radiation reaction force Fald such 
that V • V ~ 0. In a uniform constant magnetic field, the 
ALD force (7) thus reduces to 



Fald 
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where va = — 1 

(3) 

with norm va = j 
field unit vector. 

parti- 
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V is the perpendicular velocity 


( 4 ) 


B. Guiding-center transformation 

In fusion plasmas, the gyroperiod is short compared 
to the time scale associated with collisions, the ALD ra¬ 
diation reaction force, and the electric field acceleration. 
Based on this time-scale separation, the kinetic equation 
is reduced by eliminating the gyromotion in Eq. (3) us¬ 
ing Lie-transform perturbation methods [18, 19]. The 
transformation of the dissipative ALD force uses the Lie- 
transform for non-Hamiltonian dynamics, which has been 
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recently derived in Ref. [20] . In a uniform plasma, the re¬ 
sulting guiding-center distribution function for electrons 
evolves in the 2-D gyro-angle independent momentum 
space (p,^) as 


^+Vp,5-S,,4[/]=/Fp[/], (9) 

where p is the guiding-center momentum and ^ = p\\/p is 
the pitch-angle cosine. Components of the guiding-center 
momentum-space flux 


[/] = (Kfp -f Ke -f Kald) / - Dfp • Vp,if (10) 


include convective contributions from collisional drag 
Kfp, electric field acceleration Ke, and the radiation 
reaction force Kald, and a collisional diffusion tensor 
Dpp. The integral part of the guiding-center collisional 
operator is denoted /pp [/]. It describes the evolution of 
the bulk population due to collision with fast electrons. 
When momentum conservation of the electron-electron 
collision operator is essential - as for the calculation of 
electron-driven current - the term /fp[/] must be in¬ 
cluded [21]. It is generally truncated at the first order 
in Legendre expansion. The truncation ensures momen¬ 
tum conservation but allows energy dissipation such that 
it is not necessary to model energy transport to reach a 
steady-state solution. In the present paper, the term is 
set to zero unless otherwise specified. Omitting /fp[/] 
makes it possible to use the stream function to interpret 
the steady-state fluxes in momentum space, as seen in 
Sec. IV A. It also allows a comparison between the codes 
LUKE and CODE, since the integral part of the collision 
operator is not yet included in CODE. This benchmark 
is presented in Sec. IV D, where it is also shown that the 
effect of Ipp [f] on the tail of the electron distribution can 
be neglected. 

Writing out the momentum space divergence operator 
explicitly, the kinetic equation (9) becomes 


dt 


J__5 

dp 


(p^s,) 
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pd^ 




= 0 , 


( 11 ) 


where the guiding-center momentum space flux compo¬ 
nents are 


df 

Sp = —Dpp^Yp— (ATp^FP + ATp^E + Kp^ald) /, 

df 


( 12 ) 


= vT-^Dj^^Fp— -I- (ATj^e + AT^.ald) /• 


The terms contributing to these fluxes are the convec¬ 
tion and diffusion coefficients associated with the Fokker- 
Planck collision operator (which are independent of ^ for 
isotropic field particle distributions [16, 17]) 


the electric field acceleration 
Kp,E = 

and the synchrotron reaction force 

Kp^ALD =-<Jr'yp (l - , 

K\/l-^^ 

rVj^ALD =-(Jr -■ 


(16) 

(17) 
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In Eqs. (11-19), time is normalized to the collision time 
for relativistic electrons, 
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momentum p is given in units of me, the parallel elec¬ 
tric held A|| is normalized to the critical held Ec, and 
Ur = Tf-jT^r measures the relative strength (compared to 
collisional forces) of the ALD radiation reaction force 


(Ty’ 
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where Wp is the electron plasma frequency dehned by 
Wp = ^njisQm). The collisional diffusion coefficients 
^FP (p) and Rfp (p) are normalized to (mef^jTc while 
the friction coefficient Epp (p) is normalized to mc/Tc- 
An explicit form of the momentum-space huxes (12) is 
thus 


Sp = -AfP (p) ^ -f [CA|| - Fpp ip) - Grip (l - C^)] /, 




■ / Bpp (p) df 
\ P df 
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We assume cold and inhnitely massive ions, so that 
the normalized collision coefficients App, Fpp and Bpp 
only depend on Tg and [17]. Collisions with ions 
only enter the pitch-angle scattering term Bpp (p) = 
-SFP.e (p) + ^eff/(2u), while Afp(p) = d^Fpp{p)/v, 
where the normalized electron temperature is dehned as 
= ksTe/{mc^). To summarize, the normalized equa¬ 
tion (22) depends on the following independent param¬ 
eters only: the parallel electric held i5||, the normalized 
ALD frequency Gr, the electron temperature Tg, and the 
effective charge Zes- 


III. EVOLUTION OF THE ELECTRON 
DISTRIBUTION FUNCTION 


A. Force balance and runaway region 


Dpppp = App (p), 

Kppp =—Fpp {p ), 


D{{^fp = 


Bpp (p) 


(13) Some preliminary insight into runaway electron dy- 

(14) namics can be extracted from the force balance, Kp = 
Kp FP + Kp E -b Kp ald = 0, which can be expressed as 

(15) 

eA|| - Ffp (p) - Gr7P (1 - e) = 0. (23) 
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In the high velocity limit (see Appendix B) and in the 
absence of the ALD force, the force balance yields 


P 


2 


1 


(24) 


For particles with purely parallel momentum, this condi¬ 
tion determines the critical momentum pc = (i?|| — 
above which electrons are continuously accelerated. For 
particles with very large momentum p ^ 1, the condition 
(24) provides an asymptotic value that 

ATp > 0 for particles with i> ^c- 

In this paper, the runaway region is defined as the re¬ 
gion where the momentum force balance is positive, i.e. 
Kp > 0. Electrons located within this region are con¬ 
sidered runaway electrons. In the absence of the ALD 
force, this definition corresponds to the usual idea of a 
runaway electron, as the probability for an electron to 
be continuously accelerated if it enters the region where 
Kp^FP + ATp^E > 0 is very high. The momentum space de¬ 
scribed by finite-difference Fokker-Planck codes is a lim¬ 
ited domain by nature with a high-energy boundary de¬ 
fined by a maximum momentum Pmax- A proper descrip¬ 
tion of the runaway dynamics clearly requires Pmax 2> Pc- 
Then, we must distinguish between electrons located in 
the Kp^FP + Kp^F > 0 region within the code simulation 
domain p < Pmax, internal runaways, and electrons hav¬ 
ing left the simulation domain, external runaways. The 
total runaway population consists of both internal and 
external runaways, the latter being counted in the sim¬ 
ulation, albeit without following their momentum space 
characteristics. 

When the ALD force is included, and for relativistic 
electrons with p|| ^ I and p\\ ^ p±, the force balance 
(23) is approximately given by 

A|| - 1 - arpl = 0 (25) 


and yields a condition on the perpendicular momentum 
{P± = P\/l -■C^) 


P±o = 



(26) 


The momentum space in the far tail of the distribution 
is separated into the runaway region p± < p±o where 
the electric force dominates over the ALD force and the 
collisional drag such that the net force is positive, and 
a region p± > p±Q where the ALD force and collisional 
drag dominate the electric force such that the net force 
is negative. 

As the calculations in the next sections will show, in 
the presence of an ALD force, the probability for elec¬ 
trons with Kp > 0 to escape the runaway region at some 
point is high. Therefore, the concept of a runaway elec¬ 
tron in this case is more an extension of the usual def¬ 
inition than a true characteristic. We will also see that 
electrons labelled as runaways can be entirely kept within 
the simulation domain such that there are no external 
runaways. 


B. The Fokker-Planck code LUKE 

The Fokker-Plank equation is solved numerically by 
the relativistic guiding-center Fokker-Planck code LUKE. 
Equation (9) is discretized in momentum space (p, f) us¬ 
ing a 2-D finite-difference scheme with non-uniform grids 
and a 9-point differentiation procedure. A total of 1200 
points are used for the p grid, with 140 grid points de¬ 
scribing the 0 < p < 3 region with a constant grid step, 
and 1060 grid points describing the 3 < p < Pmax = 200 
region using increasing grid steps with cubic dependence. 
A total of 166 points are used for the ^ grid, with a de¬ 
creasing step size towards ^ = ±1 for increased resolution 
near the p|| axis. The code LUKE has been benchmarked 
for the usual runaway problem [11]. A benchmark in¬ 
cluding the ALD radiation reaction force is conducted 
against the Fokker-Planck solver CODE and presented 
in Sec. IV D. 

The linearization of the electron-electron collision op¬ 
erator implies that the calculation is valid only as long 
as the electron distribution is not too distorted from the 
original Maxwellian, which in practice implies that the 
runaway fraction Ur/n does not exceed a few percent. 


C. Time evolution of the distribution function 

The electron distribution evolves from an initial rela¬ 
tivistic Maxwellian distribution, which is also the steady- 
state solution of Eq. (11) for A|| = 0 and ar = 0. A 
constant electric field An = 3 is applied, the effective 
charge is Z^s = 1, the temperature is /3 = 0.1 (Tg = 5.11 
keV), and the ratio B^/n is adjusted such that ar = 0.6, 
which corresponds to typical low-density conditions in 
tokamak plasmas (i.e. n = 10^® m“^ and A = 4 T in 
the Tore-Supra tokamak). The evolution of the electron 
distribution function in the parallel direction (^ = 1) is 
shown in Fig. 1. 

In the absence of the ALD force, a runaway tail pro¬ 
gressively extends to the edge of the simulation box 
(Fig. 1(a)). The electron distribution does not converge 
to a steady-state. The runaway rate reaches an asymp¬ 
totic value (Fig. 1(e)), and the fraction of runaway elec¬ 
trons increases continuously (Fig. 1(c)). At first, the run¬ 
away rate is related to an increase in the internal runaway 
electron population. Once the runaway tail reaches the 
edge of the simulation domain (for t ^ 200), the runaway 
rate is related to the population leaving the simulation 
box and becoming external runaways. Note that, strictly 
speaking, the linearized collision operator is no longer 
valid for t/xc > 10^ as the fraction of runaway electrons 
becomes of order unity. Nevertheless, the evolution is 
continued to illustrate the absence of a steady-state so¬ 
lution. In summary, the distribution function (Fig. 1(a)) 
evolves towards an asymptotic solution with the bulk 
population depleting at a constant rate. 

In the presence of the ALD force, however, we find 
that the distribution evolves towards a steady-state solu- 
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(a) 



(b) 




(c) 


(d) 




(e) (f) 


FIG. 1. Graphs (a) and (b) represent the evolution of the electron distribution function in the parallel direction = 1) as a 
function of the electron kinetic energy graphs (c) and (d) show the fraction of runaway electrons (RE) inside and outside 
the simulation domain, and graphs (e) and (f) the corresponding runaway rates. The parameters are i?|| = 3 , ^eff = 1, /3 = 0.1, 
and Ur = 0.6. The ALD force contribution is neglected in graphs (a,c,e), whereas it is included in graphs (b,d,f). 









































6 


tion (Fig. 1(b)) as the runaway rate vanishes (Fig. 1(f)). 
No electron leaves the simulation box, and the popula¬ 
tion of internal runaways reaches an asymptotic value 
n<rln = 0.002 (Fig. 1(d)). In addition, we can observe 
the formation of a region with positive gradient in paral¬ 
lel momentum, which appears as a high-energy bump in 
the tail of the distribution function (Fig. 1(b)). Proper¬ 
ties of the time-asymptotic electron distribution function 
are examined in the next section. 



P|| 


IV. STEADY-STATE SOLUTION AND BUMP 
FORMATION 

A. Steady-state solution 

If a steady-state solution exists, it satishes the equation 
= 0. Given the axisymmetry of the momen¬ 
tum space, the divergence-free steady-state fluxes can be 
expressed as 


(a) 





27rpy^l - ^2 




(27) 


(b) 


where A(p, ^) is called the stream function. Since Sp_j • 
Vp^,cA = 0, contours of A{p,^) indicate the direction 
of the momentum-space fluxes, or streamlines. The to¬ 
tal flux of electrons between two contours is given by the 
corresponding difference in the value of A(p, ^), such that 
narrowing contours indicate regions of stronger flux. The 
stream function thus provides a very informative graph¬ 
ical representation of the steady-state fluxes in momen¬ 
tum space. While no steady-state solution to the run¬ 
away problem exists in the absence of ALD force, it is 
possible to artificially obtain a steady-state distribution 
function by adding a source term at p = 0, which com¬ 
pensates exactly for the external runaway rate. Whereas 
adding cold electrons does not change the runaway rate 
or the shape of the distribution function in the tail, it 
enables us to interpret the stream function as a represen¬ 
tation of the steady-state fluxes. 

The 2D representation of the steady-state solution 
corresponding to the simulation parameters from Sec¬ 
tion III C is presented in Fig. 2. The distribution function 
is shown in Figs. 2(a) and 2(b) for the case without and 
with ALD force, respectively. The dashed line delim¬ 
its the runaway region where Kp > 0. The correspond¬ 
ing contours of the stream function A(p, ^) are drawn in 
Figs. 2(c) and 2(d). 

In the absence of the ALD force, the runaway popula¬ 
tion peaks near the p|| axis but extends quite far in the 
perpendicular direction into the runaway region, as seen 
if Fig. 2(a). The open streamlines represented by the 
contours of the stream function in Fig. 2(c) show that 
electrons located in the runaway region are indefinitely 
accelerated and eventually escape the simulation domain. 

In the presence of the ALD force, the runaway region 
consists of a narrow band along the py axis delimited by 



(d) 


FIG. 2. Contours of the electron distribution function in 
graphs (a) and (b) and stream function in graphs (c) and 
(d) in 2-D guiding-center momentum space (p||,px) at time 
I/tc = 10®. The parameters are E|| = 3 , ZeS = 1, /3 = 0.1, 
and (Jr = 0.6. The ALD force contribution is neglected in (a) 
and (c) whereas it is accounted for in (b) and (d). 
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P± < P±o at large energies. As seen by the contours of 
the stream function in Fig 2(d), the steady-state electron 
flux is directed towards higher energies for p± < p±o and 
towards lower energies for p± > p±o. The electron tail 
is thus confined in the region close to p± = 0, as seen 
in Fig. 2(b), which creates a strong gradient in p±. As 
seen in the next section, this gradient results in strong 
pitch-angle scattering, which contributes to limiting the 
energy of runaway electrons and gives rise to a bump 
in the distribution under certain conditions. This bump 
is centered on the p\\ axis. As the minimum between 
the bulk and the bump population naturally lies in the 
vicinity of the critical field, the bump population almost 
coincides with the runaway region, such that the electron 
population can be formally separated into a bulk and a 
runaway “beam”. 


B. Perpendicular force balance 


The spherical representation (p, ^) is the natural coor¬ 
dinate system for describing collisions. Thus, it is used 
for the numerical discretization of the kinetic equation 
in the Fokker-Planck code LUKE. As seen in Fig. 2(b), 
however, the tail of the distribution function is deter¬ 
mined by the runaway region pj_ < p_Lo and the natural 
coordinate system is rather the cylindrical representation 
(P±,P||) with pj_ = pyjl - ^2 and p|| = p^. 

The transformation to the (px,P||) is detailed in Sec¬ 
tion A. Focusing on the tail region of momentum space 
near the parallel axis (characterized by p_L <g; p|| and 
P|| :§> /3), the momentum-space fluxes entering (A2) yield 
to leading order in /? 




^11 


1 + ^eff df 
2v dp± 


P\\ 


1 

-IT + (TrV 
V 


(1 + pi) 



(7rVp\ 


/, 


(28) 


where v ~ p ||/7 ~ p\\/+ p| is the velocity normalized 
to the speed of light. 

At high energy, u ~ 1 and (28) becomes approximately 


o ^ P-L fi I /"i I 2 M f 

= -r--1 -f cr^ (1 +P±j /> 

2 dp_L P|| 

•S'!! = 0-r- [pio - pi] /• 


(29) 


The parallel flux of electrons S\\ is positive for particles 
with p± < pj_o and negative for particles with p± > pj_o. 
The resulting strong perpendicular gradient enhances 
pitch-angle scattering, which creates a positive flux in 
the pj_ direction at high energy, which is illustrated by 
the streamlines in Fig. 2(d). This flux limits the exten¬ 
sion of the electron distribution to higher energies, as 
seen in Fig. 2(b). The existence of a bump in the dis¬ 
tribution is also driven by the perpendicular dynamics. 
From the expression (29) for S±, we see that the con¬ 
vective component decreases with py while pitch-angle 


scattering is independent of py for a given perpendicular 
gradient df /dp ±. Therefore, on average, electrons in the 
tail are pushed into the runaway region at lower py, while 
they are scattered away at higher py. This dynamics is 
clearly seen in the stream function plot 2(d). Naturally, 
the bump appears at the balance point where ~ 0 
and electrons accumulate as a result of the perpendicu¬ 
lar dynamics. Given the cylindrical symmetry, and since 
the perpendicular gradient is created by the parallel force 
balance, we may assume a parabolic dependence scaled 
by P_Lo around the bump location 


1 df p_L 

fdp± p^o 


(30) 


such that from (29) and for p± —>■ 0 we obtain an es¬ 
timate for the parametric dependence of the position of 
the bump 


II& oc 


1 CJr 


1 + Zeff 


iE\\ - 1) , 


(31) 


which is in agreement with the expression obtained from 
approximately solving the kinetic equation analytically 
[15]. It is quite intuitive to expect that the bump en¬ 
ergy increases with the electric field amplitude, whereas 
it decreases with the amplitude of the ALD force and the 
effective charge. However, the underlying processes are 
rather complex and involves both the parallel and per¬ 
pendicular dynamics. The parametric dependence of the 
bump location predicted by Equation (31) will be com¬ 
pared to numerical calculations in Section V. 


C. Validity of the uniform plasma approximation 


In tokamak plasmas, particles with purely parallel ve¬ 
locity are subject to an ALD force due to the toroidal and 
poloidal periodic motions. Eor a safety factor g « 1 and 
electrons with py ^ 1, the contribution from the field 
line curvature to the ALD force is derived in Appendix 
C and is expressed as (C3) 


Kr = -(Jr 



(32) 


where po = mc/{eB) is the Larmor radius of relativistic 
electrons and R is the major radius. The momentum pR 
for which the toroidal ALD and drag forces compensate 
the electric force is thus given by Uy — I — ar {po / R)^Pr = 
0, which yields 



(33) 


For the Tore-Supra example shown in Sec. Ill, po/R = 
1.5 X 10“^ and we find pr = 110, which corresponds to 55 
MeV electrons. We see in Fig. 1(b) that the combination 
of uniform plasma ALD force and pitch-angle scattering 
limits the distribution to energies much below 55 MeV. 
Toroidal effects could thus be neglected in this case. 












D. Benchmark of the solution from the LUKE and 
CODE codes 

The simulations presented in this paper were obtained 
using the code LUKE. While the code is extensively 
benchmarked for the usual runaway problem [11], LUKE 
simulations including the ALD reaction force are pre¬ 
sented for the first time in this paper. In order to bench¬ 
mark the numerical simulations, calculations from LUKE 
are compared to those from the solver CODE, which 
solves the same Fokker-Planck equation (11) but uses 
a spectral representation of the pitch-angle dependence 
[13]. The corresponding steady-state distribution func¬ 
tions are shown in Fig. 3(a) for the parameters used 
in Sec. IIIC, and two different values for the effective 
charge, Zeff = 1 and Zeff = 4. Results from the two 
codes are in excellent agreement. In particular, both 
codes show the appearance of a bump at the same en¬ 
ergy for Zeff = 1, while they show no bump formation for 
^eff = 4. 

In addition, the integral part of the collision operator 
/pp can be included in the code LUKE. Comparing the 
cases with and without /pp, we find that the distribution 
functions are very similar, as seen in Fig. 3(b). This is 
not surprising as /pp mainly affects the bulk population, 
such that it is appropriate to ignore it in the context of 
the present paper. However, as /pp ensures momentum 
conservation in the electron-electron collision operator, 
it must be included for accurate driven current calcula¬ 
tions. Indeed, the current density associated with the 
distributions shown in Fig. 3(b) is J/{ecn) = 0.015 with¬ 
out /pp while it is J/{ecn) = 0.027 when /pp is included. 
The difference arises from a shift of the electron bulk in 
the parallel direction, which is hardly visible in Fig. 3(b); 
however, the resulting asymmetry has a strong effect on 
the corresponding current. 


V. PARAMETRIC DEPENDENCES OF THE 
ELECTRON DISTRIBUTION 

The relevant physical parameters for the runaway elec¬ 
tron problem described in this paper are the electric field 
amplitude, the magnitude of the ALD radiation reaction 
force, the effective charge, and the electron temperature. 
In this section, the bump formation is characterized as a 
function of these parameters. The distribution function 
is evolved until it reaches a steady-state solution. 

In a first set of calculations, the electric field is varied 
while keeping the other relevant parameters fixed, with 
Zeff = 1, /3 = 0.1, and <Tr = 0.6. The results are shown 
in Fig. 4. We observe that the electric field must reach a 
certain threshold for the bump to appear in the tail of the 
distribution function. Above this threshold, the energy 
corresponding to the bump location increases with /1||, 
in accordance with the estimate (31). In addition we 
observe that the number of runaway electrons, i.e. the 
number of electrons with a positive parallel force balance, 



(a) 



(b) 

FIG. 3. Electron distribution function in the parallel direction 
(^ = 1) as a function of the electron kinetic energy £k- Graph 
(a) compares results from the code LUKE and GODE for 
Zeff = 1 and Zeff = 4, respectively. Graph (b) compares 
LUKE calculations in which the integral part of the collision 
operator /pp is included in the Fokker-Planck equation (9) 
to the case where it is not. Fixed relevant parameters are 
i5|| = 3, /3 = 0.1, and ar = 0.6. 


increases with the amplitude of the electric field. Note 
that the calculation was restricted to < 3.5, as the 
linearization of the collision operator fails above this limit 
since the runaway population becomes of the order of the 
bulk population. 

In a second set of calculations, the electric field is 
fixed to /?|| = 3, while we vary the amplitude of the 
synchrotron radiation force, which is proportional to B^. 
The results are shown in Fig. 5. As expected from (31), 
we observe that the bump size and the location of the 
bump maximum in energy both decrease if <Tr is in¬ 
creased, to the point where the bump disappears if ar 
is above a certain threshold. 

In a third set of calculations, the temperature is var¬ 
ied while the normalized amplitudes of the electric field 
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FIG. 4. (a) Electron distribution function in the parallel di¬ 
rection = 1) as function of the electron kinetic energy £k, 
and (b) fraction of electrons in the runaway region. The sim¬ 
ulation results are plotted for various values of the parallel 
electric field. Fixed relevant parameters are = 1, /3 = 0.1, 
and (Tr = 0.6. 
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(b) 

FIG. 5. (a) Electron distribution function in the parallel di¬ 
rection = 1) as a function of the electron kinetic energy 
£k, and (b) fraction of electrons in the runaway region. The 
simulation results are plotted for various values of the syn¬ 
chrotron reaction force amplitude. Fixed relevant parameters 
are ^eff = 1, /3 = 0.1, and £1|| = 3. 


and synchrotron radiation force are fixed. The results 
are shown in Fig. 6. We observe that the bump existence 
and energy are not affected by the electron bulk temper¬ 
ature, which is again in accordance with the analytical 
estimate (31). However, the number of electrons in the 
bump increases strongly with Tg, to the point where the 
linearization of the collision operator fails for Tg > 10 keV 
with our choice of parameters. This dependence can be 
explained by the Dreicer effect, which feeds the runaway 
population from the bulk via collisional diffusion. As 
seen in Fig. 1(b), the energy corresponding to the mini¬ 
mum between the bulk and the runaway bump decreases 
with time, until it becomes of the order of the critical 
energy. At this stage, the minimum coincides with the 
point where forces balance, such that the collisional diffu¬ 
sion in energy comes to a halt. In other words, the bump 
population increases until the negative diffusive flux as¬ 


sociated with the positive energy gradient of the bump 
is sufficient to compensate for the Dreicer flux. Since the 
latter strongly depends upon the bulk temperature, the 
bump population evolves accordingly. 

Finally, in a fourth set of calculations, the effective ion 
charge is varied while all other relevant parameters are 
fixed. The results are shown in Fig. 7. We observe that 
the bump size and energy decrease with Z^s, to the point 
where the bump disappears for > 3 with this choice 
of parameters. The effect of the ion effective charge is 
predicted by the estimate (31) and understood via the 
role of pitch-angle scattering, which is proportional to 
1 -I- Zes. For a given perpendicular gradient in the tail 
of the distribution function - which is determined by the 
parallel force balance - pitch-angle scattering is the dom¬ 
inant mechanism to extract electrons from the runaway 
region. The bump can exist only if the perpendicular 
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FIG. 6. (a) Electron distribution function in the parallel di¬ 
rection (^ = 1) as a function of the electron kinetic energy 
Sk, and (b) fraction of electrons in the runaway region. The 
simulation results are plotted for various values of the elec¬ 
tron temperature Te. Fixed relevant parameters are Z^fi = 1, 
Fl|| = 3, and Ur = 0.6. 


FIG. 7. (a) Electron distribution function in the parallel di¬ 
rection (^ = 1) as a function of the electron kinetic energy £k, 
and (b) fraction of electrons in the runaway region. The sim¬ 
ulation results are plotted for various values of the effective 
charge. Fixed relevant parameters are /3 = 0.1, E|| = 3, and 

(Tr = 0.6. 


convection due to the synchrotron radiation force and 
collisional drag dominates over pitch-angle scattering at 
the lower energies in the runaway region. 

VI. CONCLUSIONS 

In this paper, the effect of the Abraham-Lorentz-Dirac 
force in reaction to the synchrotron emission of run¬ 
away electrons is investigated for a homogeneous plasma. 
Whereas a runaway region - with positive force balance 
- can still be identified in the presence of the ALD force, 
the electron distribution decreases with momentum at 
high energy and evolves towards a steady-state solution. 
This evolution is a result of the strong pitch-angle scat¬ 
tering associated with large gradients in perpendicular 
momentum. The distribution of electrons is limited to 


energies well below the value for which the contribution 
from the toroidal parallel motion to the ALD radiation 
reaction force becomes significant in a tokamak plasma 
[9], which justifies the uniform plasma approximation to 
describe the runaway dynamics in the plasma center. 

If the electric force is large compared to the ALD force 
(proportional to B^) and the effective charge (which de¬ 
termines the rate of pitch-angle scattering), a bump cen¬ 
tered on the parallel momentum axis is formed in the 
steady-state electron distribution. It results from the 
competition between the perpendicular convection due 
to collisions and the ALD force, and pitch-angle scatter¬ 
ing. This bump encompasses almost the entire runaway 
electron population, thus formally dividing the distribu¬ 
tion into a bulk population and a runaway beam. The 
steady-state population of electrons in the bump is found 
to increase with the bulk electron temperature T^. The 
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bump size and average energy increase with the electric 
field amplitude , whereas they decrease with the ampli¬ 
tude of the ALD radiation reaction force and the effective 
charge. 

We can summarize the effect of the ALD radiation re¬ 
action force on the electron distribution in three points: 
first, in accordance with experimental observations it lim¬ 
its the energy gained by runaway electrons to the tens 
of MeV range; second, it increases the perpendicular 
anisotropy of the electron distribution, which may give 
rise to kinetic instabilities such as the EXEL or whistler 
waves [22-24]; third, it can lead to the formation of a 
bump in the electron tail, which may give rise to plasma- 
beam types of kinetic instabilities. 

Large amplitude kinetic instabilities generated by the 
runaway population could pump energy away from elec¬ 
trons and also affect their confinement. Both effects 
could be beneficial when attempting to limit the threat 
posed by runaway electrons in tokamaks. Quantifying the 
effect of kinetic instabilities requires a quasilinear treat¬ 
ment of the kinetic wave-particle interaction, which is 
beyond the scope of this paper. 

More generally, it is interesting to note that any force 
with a magnitude that increases with the particle energy 
could play a similar role as the ALD radiation reaction 
force, resulting in a maximum energy limit for runaways 
and the possible formation of a bump in the energy dis¬ 
tribution. 
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Appendix A: Cylindrical representation 
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Appendix B: High-velocity limit 

Properties of the electron distribution function are in¬ 
vestigated under the conditions that the thermal elec¬ 
tron energy is much smaller than the electron rest mass, 
namely /3 ^ 1, and that the electric field is larger than 
the critical field but much smaller than the Dreicer field 
Ejy = meaning 


1 < E|| < /3-2 (Bl) 

This ordering implies that runaway electrons are located 
in the tail of the distribution function, with a momen¬ 
tum p ^ /3 where P is the normalized thermal momen¬ 
tum. For such electrons it is appropriate to take the high 
velocity limit of the collision operator, which yields 


The kinetic equation (11) can be expressed in the 
(P||,P-l) system, which yields 


dt 


\ _^ 

P± dpA_ 


(pa5a) + ^(E||)=0, 


(Al) 


^FP {p) = (B2) 

yO 

Efp (j>) = —, (B3) 

yZ 

EFP(p) = m^. (B4) 


with the following expressions for the flux components 


Ci n n 

7-L = — Dj_|| Fp—-Dj_j_ FPT^- 

op\\ dp^ 

+ K±^Fpf + K±^Ef + A'j__ald/, 

c. n n 

^II=-^IIII.fp^-E||a,fp^ 

+ E||_fp/ -I- ATii e/ + ATii ald/, 


(A2) 


Appendix C: ALD radiation reaction force for 
purely parallel motion 

In a non-uniform magnetic field, as is found in toka¬ 
maks, the field line curvature affects the ALD radiation 
reaction force. Whereas this effect is expected to be small 
compared to the contribution from the cyclotron motion 
for particles with a significant magnetic moment, it could 
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play a role for particles with px — 0, for which the con¬ 
tribution from the cyclotron motion (8) vanishes. The 
combined effects of cyclotron motion and field curvature 
to the ALD radiation reaction force have been evaluated 
in a previous work for a purely toroidal magnetic field 
[9]. Whereas a self-consistent calculation of the ALD 
radiation reaction force in a tokamak geometry requires 
a proper guiding-center transformation [20], the impor¬ 
tance of the contribution from the field line curvature can 
be approximately evaluated by considering the motion of 
a particle with pj^ = 0. The corresponding guiding cen¬ 
ter follows the field lines with a velocity v = such 
that the ALD radiation reaction force (7) reduces to 


K = 


(Trl'^vlpl 


6-V 


{b-Vh^ +j^viib- 6 - 


(Cl) 


where the normalization of Sec. IIB is used and with 
Pq = mc/{qB). In a tokamak with major radius Rq and 
circular concentric flux-surfaces characterized by the lo¬ 


cal inverse aspect ratio e = r/i?o 1, K can be ex¬ 
pressed as 



-I-(l-b 2ecos0 {g ^ - l}) ^ -b C>(e^) , 


(C2) 


where 9 and cj) denote the unit vectors in the poloidal and 
toroidal directions, respectively, 6 is the poloidal angle, 
and q{e) is the safety factor. In the case of purely toroidal 
field lines {q oo) the results from Ref. [9] are retrieved. 

The approximation e <C 1 is valid near the plasma 
center, where in addition we typically have g ~ 1, in 
which case (C2) becomes 


K = —(7^7 



b + 0{e^) 


(C3) 
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